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Abstract. The splitting of the frequencies of p-mode multiplets 
enables information to be gained about the internal rotation of 
the sun. Such data have revealed a transition at the base of the 
convection zone from differential rotation similar to that ob- 
served at the surface to almost solid-body rotation in the ra- 
diative interior. This transition region, known as the tachocline, 
has been found to be relatively narrow and centred below the 
base of the convection zone. In this paper, the evolution of the 
transition region is investigated numerically. Without a large 
anisotropic viscosity, the depth to which it would spread in one 
solar age, under the assumption of a constant prescribed differ- 
ential rotation at the base of the convection zone, is found to be 
greater than its extent as inferred from helioseismology. In the 
second part of the paper a highly anisotropic turbulent viscos- 
ity with a large horizontal component, as suggested by Spiegel 
& Zahn (1992), is assumed. In this case, a steady tachocline is 
formed in which the advection of angular momentum balances 
the Reynolds stresses. The horizontal component of turbulent 
viscosity required to match the thickness of the tachocline to that 
obtained by helioseismology is estimated to be 5 x 1 0 5 cm 2 s" 1 . 
The transport of helium is studied in this case and is found to 
yield a sound-speed increase similar to that required by helio- 
seismology. 
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1. Introduction 

The accurate measurement of the frequencies of solar p-modes 
has enabled much information to be gained about the distribu- 
tion of angular momentum within the solar interior (Brown et 
al. 1989; Goode et al. 1 99 1 : Sekii 1991 ; Thompson et al. 1996). 
This has made the search for a good theoretical understand- 
ing of the processes at w ork within rotating stars all the more 
important. 

The possible existence of a meridional circulation, w hich 
could transport material and angular momentum within the so- 
lar radiative interior, was first proposed by Eddington in the 
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early part of this century. However, a consistent picture of the 
distribution of angular momentum to w hich the sun evolves has 
been difficult to ascertain. Tassoul & Tassoul (1982), amongst 
others, argued that the sun must evolve to a state where a balance 
exists between the transport of angular momentum by advection 
and its transport by turbulent viscosity. This turbulence could 
very well arise from shear instability (Spiegel & Zahn 1970), 
with other instabilities almost certainly playing a role. 

A particularly interesting aspect of the solar angular- 
momentum distribution is the transition layer between the 
rigidly-rotating radiative interior, and the differentially-rotating 
convection zone. This layer, known as the tachocline, is be- 
lieved to be an important component of the dynamo which is 
responsible for the solar cycle - its strong shear has the ability 
to transform a poloidal field into a strong toroidal field, which 
can rise up buoyantly through the convection zone and erupt on 
the solar surface as sunspots (Spiegel & Weiss 1980). 

In this paper, the time evolution of differential rotation in 
the tachocline is studied by means of a numerical simulation, 
assuming the rotation rate at the base of the convection zone 
to be prescribed. The equations solved are similar to those of 
Spiegel & Zahn (1992). Various other simplifying assumptions 
are made, the details of which are described in Sect. 2.2. It is 
initially assumed that the flow is not highly anisotropic - in this 
case the differential rotation is found to spread far into the radia- 
tive interior. A highly anisotropic viscosity is then introduced in 
Sect. 5 - this has the effect of limiting the spreading of the dif- 
ferential rotation and establishing a balance between advection 
of angular momentum and transport by turbulent viscosity. 

The effect of the corresponding circulation on the transport 
of helium is studied in Sect. 6. Helium originating in the convec- 
tion zone, which has settled under gravity into the tachocline, is 
mixed back into the convection zone. This causes a reduction in 
the mean molecular weight in the tachocline, and a consequent 
increase in the sound speed there. Since helioseismic inversions 
in this region indicate a higher sound speed in the sun than in 
current models, mixing improves the agreement with such in- 
versions. 
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2. Governing equations 

2.7. Equations in a rotating frame 

In order to describe differential rotation in the solar interior it is 
useful to write the equations of motion with respect to a frame 
rotating with angular velocity L2. The equations of conservation 
of mass, momentum, and energy become 


dp 

dt 


+ V • (pv) = 0 , 


dv 

— + (v ■ V)r + 21 1 x v 
dt 


(1) 


= - Vp - pV<L + V ■ | |r 1 1 , (2) 


pT^-+v-VSj = V ( X VD , (3) 

where p, T and p are the density, temperature and pressure 
respectively, v is the fluid velocity, S is the specific entropy, \ 
is the thermal conductivity, is the gravitational potential, and 
||r|| is the turbulent stress tensor. 


2.2. Assumptions 

The equations are solved under the assumption that the mean 
structure and the large-scale flow are axisymmetric about the 
axis of rotation (the small-scale turbulence is certainly not). The 
system can therefore be described solely in terms of the spherical 
coordinates r and <9, where 9 is the colatitude. The solution is 
sought in a spherical shell with the outer edge corresponding to 
the base of the convection zone. 

All structure quantities, for example the pressure, p, are sep- 
arated into a mean value on the sphere plus a perturbation: 

p(r, 9 . t) = p(r, t) + p(r, 9 , t) , (4) 

where f p(r 7 9. t) sin 6 d9 = 0. The angular rotation rate, lc, is 
written as 

u ;(r, 9 , t) = ■ - l ~- = u;(r. t) + u)(r. 0, 7) , (5) 

r sin 0 

w r here / T’(r, 0, i ) sin 3 9d9 = 0. 

Equations ( 1 )-(3) are solved under assumptions which are 
very similar to those employed by Spiegel & Zahn (1992). They 
are as follows: 

1 . The time scale of the flow is very much longer than the sound 
crossing time, and the circulation can therefore be calculated 
using the anelastic approximation, i.e. that V-(pv) = 0. This 
cannot hold strictly as it would imply dp/dt = 0 which is 
clearly not the case. However, since the time scale of the 
circulation is very much shorter than the time scale on w hich 
the density is modified (the Eddington-Sweet time scale). 
dp/dt may safely be neglected. 

2. Since the time scale of the flow is long compared to the ro- 
tation time scale, the inertial and viscous forces are small 
compared to the Coriolis forces in the radial and latitudi- 
nal directions; in these two directions a geostrophic balance 


exists between the Coriolis force and the pressure gradient. 
Since no zonal pressure gradient can exist owing to the as- 
sumption of axisymmetry, the azimuthal Coriolis force must 
be balanced by inertial and viscous forces. 

3. The oblateness in the figure of the sun caused by the cen- 
trifugal force is small enough that the level surfaces of the 
effective potential can be taken to be spheres, with a com- 
pensating thermal source term added to the energy equation 
(Zahn 1992), 

£es = ^jEqP 2 (cos9) , (6) 

M 

w'here L and M are the solar luminosity and mass respec- 
tively, Pi denotes the second Legendre polynomial, and in 
the case of uniform rotation, Eq may be represented as 


Eq = 2 




9 

9 


(7) 


where Q = |fi|. In the outer regions of a star, the non- 
spherically symmetric component of the gravitational ac- 
celeration may be represented as 



4. Second derivatives of thermodynamic quantities are ne- 
glected, i.e. quantities like the specific heat at constant pres- 
sure, C p , and the thermal conductivity, are assumed to be 
constant. This leads to a linear, ideal-gas equation of state 


V _ T + p 
P T p' 


(9) 


No explicit assumptions are made regarding the thickness, 
h, of the shell in which the simulation is performed. However the 
accuracy of the simulation is sensitive to h through the neglect 
of modifications to the background state caused by the changing 
1 2(r, t). Since the modification of the background state involves 
a displacement of the order of hZJ/Q, w hich in turn reacts back 
on2(r, t) by an amount hZJ/r , this approximation becomes pro- 
gressively worse as the thickness of the shell increases. If the 
modification of the background state were included, the effect 
would be to reduce the pressure fluctuation, p, and thereby also 
to reduce the predicted circulation velocity. 


23 . Simplified equations 

Under the above assumptions, the radial component of the mo- 
mentum equation reduces to 

2r sin 2 9 HT' =- , (10) 

p or p 

where g is the acceleration due to gravity, with g assumed to 
be a constant, prescribed function of radius, since the layers 
of interest are at a large distance from the solar core. The 9 
component of the momentum equation becomes 


2 r sin 9 cos 9 


1 dp 
rfid9 * 


( 11 ) 


f 
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The o component of the momentum equation, which expresses 
the balance between Coriolis, inertial, and viscous forces, be- 
comes 


. _ 3va dVft, „ 

r sin 9 ~~~ + i‘ r — + v$ + 2fi (ty sin 9 + v$ cos 6) = 
at or du 


sin 2 9 d / 4 N 

pr- or \ or 


1 


9 


. 3 

rsi„=ea9\‘ ,,s "' < 'a? 


( 12 ) 


where and t/y are the horizontal and vertical components of 
the turbulent viscosity. 

As a result of the anelastic approximation, the circulation in 
the spherical shell may be represented by a streamfunction 


assumed to be prescribed by the convection zone above, and to 
have a form corresponding to that inferred from helioseismol- 
ogy. An inversion has recently been carried out for the solar 
internal rotation using splitting data from the GONG network 
(Thompson et al. 1996). At the base of the convection zone, this 
inversion leads to the following expression for £2 bcz : 

— ^ = 456 - 72 cos 2 9 - 42 cos 4 9 nHz . (16) 

This is used to give the outer boundary condition on u 

Q + jj = fi bcz at r = r bcz , (17) 


_ 2 ■ * d * 

pr a»«*v= gj. 

_ ■ , d * 
pr sin 9v0 = - — , 
or 


(13) 

(14) 


where ty a; : e are, respectively, the radial and tangential com- 
ponents of the velocity. Since the system is assumed to be ax- 
isymmetric about the axis of rotation as well as having reflective 
symmetry with respect to the equatorial plane (a consequence of 
the symmetry of the imposed boundary conditions), v$ should 
vanish at 9 = 0 and 9 = tt/2. Eq. (14) implies that 'P must 
also vanish at these values of 9. Integrating Eq. (13) with re- 
spect to 9 , this condition corresponds to the requirement that 
f ty sin 9d9 = 0, or that there is no net inward or outward mass 
flux. 

The energy equation consists of a spherically-averaged part, 
neglected in this paper, and a perturbation given by 


dT 

~dt 


-= d In p 1 

-tyr— AVad-V)+ — 

dr P C P 


’ 1 d 

( 

r 2 dr 



1 


r 2 sin 9 d9 


n dT 


+ P^ES 


(15) 


where r bc2 is the radius of the base of the convection zone. Q is 
chosen to be the rotation rate of the radiative interior below the 
tachocline - helioseismic inversions indicate that this value is 
equal to the rotation rate of the convection zone at a latitude of 
about 30°, corresponding to Q/2tt % 437nHz. 

The second outer boundary condition is derived from the 
requirement that the partial derivative of the rotation rate with 
respect to radius be continuous. Since helioseismic inversions 
indicate that the rotation rate depends little on radius deep within 
the convection zone, this boundary condition becomes: 

7^=0 at r = r ba . (18) 

The last two boundary conditions, 

O 9 ) 

are imposed at the inner edge of the region of solution. 


3. The solution 

3.1. numerical method 


where V is the temperature gradient d In T/d In p, and V ad is the 
thermodynamic derivative 3 In T/d In p at constant specific en- 
tropy. The potential singularity in ty caused by V ad - V passing 
through zero is avoided by locating the base of the convection 
zone at the point where the convective heat flux becomes negli- 
gible. Owing to convective penetration, this is below' the point 
where V = V ad , and the stratification is therefore subadiabatic 
in the w r hole computational domain. 

After about one thermal relaxation time, a balance is struck 
between the terms on the right-hand side of this equation, and 
the time derivative on the left-hand side may be neglected. It is 
not strictly zero, but since the system evolves on a time scale 
comparable with the local Eddington-Sweet time (Spiegel & 
Zahn 1 992). it is small enough that the terms on the right balance 
almost perfectly. 

2.4 . Boundary conditions 

Boundary conditions are imposed at the inner and outer edges 
of the region of solution. The system is fourth order and there- 
fore four boundary conditions are required. The rotation rate is 


The evolution of the rotation rate u ;(r. 9 , t) with time is studied 
numerically. Eqs. ( 10)-( 15) are solved using finite differences to 
represent partial derivatives with respect to r and 9 , and using a 
first-order explicit scheme to carry out the evolution in time. A 
uniformly spaced grid, with spacings Ar and A 9, is set up with 
equal numbers of points in the r and 9 directions, r is chosen to 
vary Tom r bcz /2 to r bcz , while 9 varies from 0 to tt/2. 

Ir itially. Eq. ( 1 1 ) is integrated to find p(r, 9. t) in terms of 
^’(r, t . t). Eq. (10) is then solved to find p(r, 9 . t). T(r , 9. t ) ma\ 
then 1 e obtained from the equation of state. After neglecting the 
time derivative in Eq. (15), the radial component of the velocity, 
ty rmy be found by taking the Laplacian of T(r. 9 , t). 

Tie streamfunction ^(r. 9. t) is now evaluated. This is done 
by int ^grating Eq. ( 1 3) with respect to 9. The tangential compo- 
nent ( f the velocity, may then be evaluated using Eq. (14). 
Finally, the rotation rate is updated using Eq. (12). The advec- 
tive Urms (the second and third terms on the left-hand side of 
this equation) are included, but turn out to be small in compari- 
son with the Coriolis term (the last term on the left-hand side), 
since he Rossby number (the ratio of the advective term to the 
Corio is acceleration) is small (Spiegel & Zahn 1992). 


f 
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cj (r,6) ^(r t 0) u(r,9) ^(r,6) 




Fig. 1. The radiative spreading of differential rotation imposed at the 
edge of the radiative zone after one solar age. The left-hand quadrant 
shows contours of w(r, 0) with solid contours denoting prograde dif- 
ferential rotation, while the right-hand quadrant shows contours of the 
streamfunction ^(r, 0), with solid contours corresponding to clock- 
wise circulation. The rotation of the base of the convection zone is 
here chosen to match that of the deep radiative interior at a latitude of 
approximately 30°. 


The stability of this explicit time-evolution scheme is deter- 
mined by the condition that the time step, A f, should be shorter 
than the time for the differential rotation to spread radially a 
distance Ar, 


At<T Es ^j , 

where r ES is the local Eddington-Sweet time. 


( 20 ) 


3.2. Reference model 

The spherically-averaged values of p, p and X, along with the 
values of C p and g , are derived from a reference solar model 
(7 is assumed to take the ideal-gas value of 5/3). This reference 
model is constructed using the MHD equation of state (Mihalas 
et al. 1988) and the OPAL opacities (Iglesias & Rogers 1991) 
with an assumed heavy-element abundance of Z = 0.02. The 
assumed solar age is 4.6 x 10 9 years. Since the outcome of the 
current calculation is not particularly sensitive to the details of 
the reference model, discussion of other details of its construc- 
tion are deferred. 


Fig. 2. As in Fig. 1, but with the boundary condition on u; at the base 
of the convection zone given in this case by Eq. (21). 


of angular momentum, the last term in Eq. (12), may be ne- 
glected in comparison with the vertical transport as long as the 
tachocline is not too thick. Assuming the vertical transport term 
to be small also, the result of the numerical calculation described 
is shown in Fig. 1. 

As can be seen, after one solar age, the transition region 
from differential rotation to rigid rotation has spread signifi- 
cantly into the radiative interior. Another feature of note is that 
the contours of the streamfunction 'T, shown in the right-hand 
quadrant, are nearly parallel to the axis of rotation at the edge 
of the region of solution. This is because the imposed rotation 
law at the base of the convection zone forces the streamlines 
to be parallel to lines of constant specific angular momentum; 
since the latter are almost parallel to the axis of rotation, owing 
to the ratio AQ/Q, being relatively small, the streamlines must 
also be nearly parallel to the axis of rotation. 

A second calculation is performed under the assumptions of 
the work of Spiegel & Zahn (1992); in this analytical study, the 
radial component of the velocity was neglected in determining 
the evolution of the differential rotation, on the grounds that in 
a thin shell it would be small in comparison with the tangential 
component. This approximation allowed the equations to be 
separated in the r and 0 directions, with the caveat that uj was 
forced to be zero at 9 = rr/2. To test the effect of this restriction, 
a numerical simulation is carried out choosing 


jj ( r Ecz . 0 . t ) = —72 cos 2 0 ~ 42 cos 4 0 nHz . (21 ) 


3.3. Initial conditions 

The initial conditions chosen for all simulations are *c(r. O.t - 
0) = 0. Since there is no w ay of knowing what the actual initial 
conditions should be in the solar case, there is an implicit as- 
sumption that the results obtained are not particularly sensitive 
to the initial conditions. 

4. Radiative spreading 

In this section, the turbulent viscosity is assumed not to be highly 
anisotropic. As a consequence of this, the horizontal transport 


in accord w ith Eq. (16). The result of this second simulation is 
shown in Fig. 2. 

The depth to which differential rotation imposed at the edge 
of the radiative interior would diffuse in one solar age is seen 
to be roughly 0.2 R t in the cases shown in Fig. 1 and Fig. 2. 
This is somewhat smaller than the value obtained by the an- 
alytical calculation of Spiegel & Zahn (1992), namely about 
0.26257? 3 . but still conflicts with the helioseismic determina- 
tion of tachocline thickness and depth carried out by Koso- 
vichev (1996). In the latter paper, a depth of 0.692 ± 0.0057? 7 
and a thickness of 0.09 ± 0.047? r were obtained by fitting a 



1226 


J.R. Elliott: Aspects of the s< »lar tachocline 



Fig. 3. The function A 3 (r) for the cases shown in Fig. 1 (dashed line) 
and Fig. 2 (dot-dashed line), and the form deduced by helioseismology 
(solid line). 


differential-rotation profile parameterized by depth and thick- 
ness to the solar p-mode splitting coefficients. In order to con- 
trast the results of the radiative spreading described here with 
the results of helioseismology, the function A 3 (r), introduced 
by Kosovichev (1996) as the coefficient of the associated Leg- 
endre polynomial P^icosO) in an expansion for the rotation 
rate, is calculated for the two cases described in this section 
and compared with his parameterized best fit to the solar p- 
mode splittings. This particular coefficient was investigated by 
Kosovichev (1996) on the grounds that it is especially sensitive 
to the variation of angular rotation rate in the solar tachocline. 
Fig. 3 represents this comparison; the dashed line shows the 
case of interior rotation matching convection-zone rotation at a 
latitude of 30°, the dot-dashed line shows the case where u is 
zero at the equator, while the solid line shows the best fit to the 
helioseismology data determined by Kosovichev (1996). 

The spreading predicted by the results of these simple cal- 
culations is seen to be much greater than that determined ob- 
servationally. In addition, it was here assumed that Q, was con- 
stant over the whole life time of the sun. In reality, the sun has 
been slowing down due to the braking effect of the solar wind, 
and was rotating faster in the past. The spreading would be in- 
creased if this slowing down was taken into account, making the 
discrepancy between the predictions of this model and obser- 
vations even greater. Some way needs to be found of inhibiting 
the spreading of the tachocline into the radiative interior. 


5. Anisotropic turbulent viscosity 

Spiegel & Zahn (1992) suggested that strongly anisotropic tur- 
bulence could inhibit the spreading of the transition region into 
the radiative interior of the sun, and that in the solar tachocline 
a balance may exist between the advection of angular momen- 
tum by the circulation and the Reynolds stresses acting on the 
horizontal shear. Since the thickness of the tachocline would 
depend on the horizontal component of the turbulent viscosity, 
a knowledge of the thickness would enable this component to 
be estimated. 





Fig. 4. The steady state reached when the advection of angular mo- 
mentum by the meridional circulation balances the Reynolds stresses. 
The horizontal component of the turbulent viscosity is here chosen to 
be 5 x 10 4 cm 2 s _I , while the ratio vh / ia / = 1000. 


When 



( 22 ) 


as may be the case when the stable stratification of the subcon- 
vecti\ e layers produces strongly anisotropic turbulence, the first 
term on the right-hand side of Eq. ( 1 2) may be neglected to give: 


• / c) Ufa d ijsi) 

r sin 6 — + t’ r — — + v$ - = -2f> (v r sin 0 + v e cos Q) 
at dr du 


1 


a 


1 2 QdO 


e t0 ■ <23 > 


Once a balance has been established between the Coriolis 
and v scous forces (the two terms on the right-hand side of 
this equation), a steady state is reached. An example of such a 
stead > state is shown in Fig. 4, which corresponds to a horizontal 
turbulent viscosity of 5 x 10 4 cm 2 s -1 , and a ratio j/hA'v = 1000. 
This ratio is sufficiently large that the criterion given by Eq. (22) 
is met 

Tf e rotation rate at the base of the convection zone matches 
that oJ the deep radiative interior at a latitude of approximately 
40°, much as predicted by the analytical calculation of Spiegel 
and Z^hn (1992). The circulation also shows the octopolar con- 
figuration that they described, with two cells in each hemisphere. 
The thickness of the tachocline may be compared with the ap- 
proximate formula given in that paper. 


h = 2 C 000 (k/V h ) 4 • (24) 

where k = \/pC p . Given - 5 x I0 4 cm 2 s _1 and using the 
value of k from the base of the convection zone, this formula 
predic s a thickness of about 90 000km, or 0.1 3R- . This is 
somew hat larger than the tachocline thickness shown by our nu- 
merical calculation for this value of (approximately 0. 1 R -,), 
The discrepancy reflects the fact that the approximate formula 
(24) was obtained under the assumption that the tachocline is 
vanish ngly thin and that a/ varies much more rapidly with ra- 
dius th in do other structure variables; when allowance is made 


l 
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Fig. 5. The function A 3 (r) in the case of strongly anisotropic turbulent 
viscosity. The dashed line corresponds to ^ = 5 x 10 4 cm 2 s _1 , the 
dot-dashed line corresponds to uh = 5 x 10 5 cm 2 s _1 , while the solid 
line corresponds to the form deduced by helioseismology. 

for the variation of the background state with depth, as in this 
numerical calculation, the tachocline is found to be thinner than 
predicted by Eq. (24). 

Fig. 5 shows the comparison between the function ,4 3 (r) 
for two models with strongly anisotropic turbulent viscosity 
(but different values of this viscosity) and that derived by Koso- 
vichev (1996) from helioseismology. The dashed line shows 
Ai(r) corresponding to the case shown in Fig. 4, while the dot- 
dashed line shows the same function for a horizontal turbulent 
viscosity ten times larger (5 x 1 0 5 cmV 1 * ). The solid 1 me shows 
the value obtained from helioseismology. 

The better fit to the helioseismically deduced A 3 (r) is ob- 
tained with the higher value of 5 x 10 5 cm 2 s _1 for the horizontal 
turbulent viscosity plotted in Fig. 5. This value is significantly 
smaller than the lower bound of 1.5 x 10 8 cm 2 s“ 1 stated by 
Zahn (1992) based on Eq. (24). If a tachocline thickness of 
0.1/20 or 70000km were presumed, then Eq. (24) would re- 
quire the ratio k / to be approximately 150. n is approximately 
2 x 10'crrrs ! at the base of the convection zone, leading to 
i/H ^ 10 5 cm 2 s“ 1 . The estimate of 1.5 x 10 8 cm 2 s“ l as a low'er 
limit lor i/ H would thus appear to be wrong, as it would lead 
to a tachocline several times narrower than the resolution of 
helioseismic inversions for the rotation rate. 

Zahn (1992 ) used this limit to argue that if all the energy go- 
ing into turbulent motions (as calculated from the shear implied 
by helioseismic measurements) were dissipated by viscosity on 
a small scale, then the vertical component of the turbulent vis- 
cosity required would imply a turbulent diffusivitv too large to 
tolerate the observed surface abundance of lithium. This find- 
ing may simply be a consequence of the overestimate of the 
horizontal component of the turbulent viscosity. 

The thickness of the tachocline in Fig. 4 shows a clear vari- 
ation with latitude. There is some evidence in helioseismic in- 
versions for rotation rate (Thompson et al. 1996) for such a 
variation being present in the sun, giving additional support to 
this theoretical description of the solar tachocline. The actual 
magnitude of the variation has yet to be quantified, and awaits 
more accurate determination of p-mode frequency splittings. 


Under the assumption of no net torque at the base of the 
convection zone, models of the solar tachocline with strongly 
anisotropic viscosity predicta value of about 41 6nHz % 0.9 lft e 
for the rotation rate of the solar radiative interior (f> e is the 
equatorial rotation rate); the value is determined solely by the 
imposed rotation law at the base of the convection zone, a de- 
pendence which has been described analytically by Spiegel & 
Zahn ( 1 992). This differs from the current best helioseismic esti- 
mates, e.g. from GONG data, of about 435nHz ^ 0.95f2 e . Thus, 
in this model, the rotation rate of the radiative interior matches 
that of the convection zone at a latitude of about 40°, while he- 
lioseismology indicates that the two are equal at a latitude of 
nearer 30° in the sun. 

In the presence of a torque, such as could be produced by 
angular-momentum loss associated with the solar wind, it might 
be expected that a large differential rotation would be set up, as 
suggested by Zahn (1992). Since helioseismic inversions indi- 
cate that there is no such differential rotation in the radiative 
interior, some mechanism must inhibit this process. It has been 
suggested by Rosner & Weiss ( 1 985 ) and Mestel & Weiss ( 1 987) 
that a weak connected magnetic field may cause the solar interior 
to rotate uniformly. Alternatively, Zahn (1994) has suggested 
that gravity waves excited by penetrative convection may trans- 
port angular momentum radially and help to enforce uniform 
rotation. The relative importance of such processes may become 
clearer as more detailed helioseismic data become available. 

6. Mixing in the turbulent tachocline 

Recent helioseismic inversions for the sound speed in the 
tachocline show a higher (Sc 2 fc 2 % 0.004) value in the sun 
than in the latest models. These models incorporate gravita- 
tional helium settling (Michaud & Proffitt 1993), which causes 
the helium abundance to increase with depth below the base of 
the convection zone. It has been suggested that the discrepancy 
could be due to the circulation in the turbulent tachocline mixing 
helium back into the convection zone, thus reducing the mean 
molecular weight in the tachocline and thereby increasing the 
sound speed. 

In order to test this theory, the evolution of helium abundance 
is investigated using a model incorporating both gravitational 
settling and advection of helium. The circulation responsible for 
this advection is assumed to be that associated with a steady- 
state tachocline having = 5 x 10 > cm 2 s” ! ; it is additionally 
assumed that the circulation has not changed w ith time. A diffu- 
sion term is included to make the solution smooth at the interface 
between the tachocline and the convection zone: the diffusion 
coefficient drops exponentially with depth below the base of 
the convection zone (with e-folding distance 0.01/2 r ). where 
it has a value of 200cm : s~ 1 . The horizontal turbulent diffusion 
{corresponding to u H ) of helium is not included. The mass of the 
convection zone is assumed to decrease linearly with time up to 
the present-day solar age from an initial value of 0.02564/ : to 
a present-day value of 0.01834 /t . This linear trend, as well as 
the initial and final masses, are derived from the evolution of a 
standard solar model. 
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whi :h depends in turn on the chemical composition (a fully- 
ioni zed ideal gas is assumed). After convolution with Gaussian 
func tions chosen to have widths reflecting those of optimally- 
localized averaging kernels used in helioseismic inversions, the 
dashed curve is obtained. The large hump in Scr/c 2 below the 
base of the convection zone is very similar to that observed in 
helioseismic inversions for sound speed, and has an approxi- 
mately correct amplitude (0.004). The mixing of helium back 
into the convection zone thus appears to be a good candidate 
for explaining the sound speed discrepancy in the region of the 
tachocline. 


Fig. 6. The spherically-averaged hydrogen profiles obtained without 
advection of helium (solid line) and with advection (dashed line). 



Fig. 7. The difference in sound speed between a model with advection 
by the tachocline circulation and a model without such advection (solid 
line); the diamonds {joined by the dashed line) show this quantity con- 
volved with Gaussian functions having widths corresponding to those 
of optimally-localized averaging kernels from helioseismic inversions. 

The result of the calculation is presented in Fig. 6. The solid 
line shows the spherically-averaged hydrogen abundance, A' 
(after subtraction of the initial abundance, A"o), without advec- 
tion by the circulation, while the dashed line shows the same 
quantity when advection is included. Without advection, the in- 
crease in convection-zone hydrogen abundance by the present 
day, 0.032, agrees reasonably well with that given by more de- 
tailed evolution calculations that include helium settling. When 
advection is included, helium is mixed back into the convection 
zone and the increase in convection-zone hydrogen abundance is 
reduced to 0.024. Advection also removes the steep gradient of 
hydrogen abundance just below the base of the convection zone, 
thus increasing the sound speed in this region and improving the 
agreement with helioseismic inversions. Looking more closely 
at the shape of the dashed line, the hydrogen abundance is seen 
to descend from the convection-zone value via two plateaus, 
which correspond to the two circulation cells seen in a radial 
cut through the streamfunction shown in Fig. 4. 

The solid line in Fig. 7 shows the quantity 8r 2 jc 2 , reflecting 
the difference in squared sound speed between the model w ith 
advection and that without. This quantity is assumed to be pro- 
portional to the fractional difference in mean molecular weight. 


7. Conclusion 

This paper has confirmed some of the findings of Spiegel & Zahn 
(1992). Without anisotropic viscosity, the differential rotation 
of the convection zone would diffuse into the radiative zone at 
too great a rate to comply with observations. When anisotropic 
viscosity is introduced, the spread of the differential rotation 
is inhibited, and a steady-state is reached in which advection 
of angular momentum balances Reynolds stresses. The rotation 
rate at the base of the transition region so formed is equal to the 
rotation rate of the base of the convection zone at a latitude of 
approximately 40°. 

The simple viscous stress operator considered here is not 
the only case one could imagine. For instance, the coefficients 
of tu-bulent viscosity may vary with latitude, since turbulent 
motions are more strongly rotationally constrained at the poles 
than at the equator. 

While the results presented here are mostly in agreement 
with ihose of Spiegel & Zahn ( 1 992), the estimated value of the 
horizontal component of the turbulent viscosity, 5 x 10 5 cm 2 s“' , 
is significantly smaller than the lower bound given by Zahn 
( 1 992), and is only a few- orders of magnitude larger than the mi- 
crosc opic viscosity. This implies very short characteristic length 
scale i ( 10 4 cm) at the velocities associated w ith convective mo- 
tions. or very low velocities (10~ 5 cms“ l ) at the length scale 
associated with the tachocline. It should be noted that given the 
uncertainty associated with current estimates of the tachocline 
thickness, and the weak dependence of this thickness on the hor- 
izontal turbulent viscosity, uu could in reality be significantly 
largei than the value obtained here. 

S nee the circulation associated with the steady-state 
tachoriine advects helium, a certain amount of mixing occurs 
below the base of the convection zone. Numerical calculations 
taking into account the interaction between this mixing and 
graviiational settling show- that the sound speed could be mod- 
ified in such a way as to remove the discrepancy found in he- 
liosei mic inversions. This gives independent confirmation of 
the thickness of the tachocline as deduced from rotation-rate 
inversions. 

In obtaining these results, several overall simplifying as- 
sumptions have been made. It has been assumed that the con- 
vection zone sets the boundary conditions for the radiative in- 
terior, which does not allow for any feedback of the radiative 
interic r on the convection zone. The effects of convective over- 
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shoot have not been considered either - these would affect the 
boundary between the radiative interior and the convection zone 
over a region comparable with the pressure scale height. Tur- 
bulent transport of heat has also been neglected in comparison 
with radiative transport, an approximation which is valid so long 
as the tachocline is not too thin. A better model would consider 
the full interaction between the convection zone and radiative 
interior, allowing for the effects of penetrative convection and 
turbulent transport of heat - such a model w'ould certainly be 
much more complex than that considered here. 

What has emerged from the results presented here is the 
diagnostic potential of rotation-rate inversions using frequen- 
cies obtained from such experiments as MDI and GONG. By 
measuring the thickness of the tachocline, the variation of this 
thickness with latitude, and the interior rotation rate, many as- 
pects of the dynamics of this fascinating region of the sun may 
be studied. The great improvement which may be expected in 
the accuracy of helioseismic measurements over the next several 
years will enable a much better understanding of the tachocline 
to be reached. 
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